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ABSTRACT 


This  report  presents  one  approach  for  modeling  smallpox  and  inhalation  an¬ 
thrax  outbreaks  using  ODEs  (ordinary  differential  equations).  This  approach  is 
related  to  a  standard  SEIR  (susceptible  exposed  infected  recovered)  model.  For 
each  model,  we  define  the  states  that  characterize  the  uninfected  and  infected  pop¬ 
ulations,  the  parameters  governing  disease  progression,  and  the  ODEs  that  govern 
the  transitions  between  the  population  states.  In  both  models,  medical  capacity 
and  treatment  limitations  are  considered.  To  quantify  the  benefit  of  an  early  pub¬ 
lic  health  response,  the  number  of  cases  and  deaths  resulting  from  an  outbreak  are 
determined  as  a  function  of  delay  in  public  health  response.  The  smallpox  model  in¬ 
dicates  that  early  initiation  of  a  mass  vaccination  campaign  can  significantly  reduce 
the  number  of  deaths.  The  anthrax  model  indicates  that  distribution  of  antibiotics 
at  a  high  rate  within  the  first  day  following  a  large  attack  can  save  nearly  all  those 
exposed.  Future  work  will  focus  on  replacing  the  ODEs  with  probability  distribution 
functions  based  on  data  from  outbreaks;  doing  so  will  lead  to  a  more  accurate  model 
of  the  incubation  periods  and,  in  turn,  a  more  accurate  estimate  of  the  benefit  of 
an  early  response. 
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1.  INTRODUCTION 


This  report  presents  one  approach  for  modeling  smallpox  and  inhalation  anthrax  outbreaks 
using  ODEs  (ordinary  differential  equations).  This  approach  is  related  to  a  standard  SEIR  (sus¬ 
ceptible  exposed  infected  recovered)  model;  see  [2],  for  example.  One  benefit  of  such  a  simple 
framework  is  greater  transparency  and  fewer  model  parameters  that  must  be  estimated  and  varied. 
Other  types  of  models  include  more  complicated  ODE  formulations  [35]  and  stochastic  formulations 
with  homogeneous  mixing  [20]  and  heterogeneous  mixing  [6,11-];  homogeneous  mixing  assumes  all 
people  in  the  susceptible  population  are  equally  likely  to  become  infected.  See  [8,24]  for  a  general 
review  of  these  approaches. 

The  two  ODE  models  are  based  on  a  framework  of  states  described  in  section  2.  For  each 
model,  we  define  subcategories  of  the  primary  states,  the  disease  progression  parameters,  and  the 
transition  rates  between  the  states.  In  both  models,  medical  capacity  and  treatment  limitations 
are  considered.  Smallpox  is  modeled  as  a  single-stage  illness  with  capacities  defined  for  quarantine. 
Inhalation  anthrax  is  modeled  as  a  three-stage  illness  to  which  different  medical  treatments  and 
Ccipacities  are  applied. 
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2.  DISEASE  MODEL  FRAMEWORK 


The  populations  of  both  models  can  be  separated  into  the  six  primary  states  shown  in  Table  1. 
For  each  model,  these  states  are  further  broken  down  into  additional  states. 


TABLE  1 

Basic  states  of  the  disease  models,  which  are  functions  of  time 


Disease  States 

Notation 

Unexposed  at  time  t 

U{t) 

Exposed  people  at  time  t 

Eif) 

Symptomatic  people  at  time  t 

S{t) 

Vaccinated  people  at  time  t 

Vit) 

Recovered  people  or  people  with  immunity  at  time  t 

Rit) 

Dead  at  time  t 

D{t) 

Transitions  between  these  states  are  governed  by  ODEs.  Let  P{t)  E  be  a  vector  of  the 
number  of  people  in  each  state  at  time  t,  so  that  Pk{t)  is  the  number  of  people  in  the  fcth  state 
of  the  model  at  time  L  Given  the  initial  conditions  P(^o)  ^md  the  differential  equations,  we  can 
use  Euler’s  method  to  estimate  the  number  of  people  in  each  state  at  time  Euler’s  method  is 
a  typical  method  for  solving  differential  equations,  using  derivative  information  from  the  previous 
time  step  Ut-i  to  estimate  the  number  of  people  in  the  kth  state  at  time  tn- 

P/ci^n)  —  Pki'^n-l)  +  F^(^71-i)AL  (1) 

where  Pl^{t,t-[)  is  the  derivative  of  the  kth  population  at  time  t^-i.  For  convenience,  we  define 

=  (2) 

where  is  the  rate  of  people  moving  from  the  ith  into  the  kth  state  at  time  t,  X  is  the  set 

of  states  that  flow  into  the  fcth  state,  and  J  is  the  set  of  states  into  which  the  kth  state  flows 
{XnJ  =  0).  When  the  ?'th  state  does  not  flow  into  the  jth  state,  —  0  for  all  t.  For  notational 

convenience  when  defining  the  rates  we  replace  the  notation  with  and  t  and  j 

with  the  notation  for  the  states  shown  in  Table  2. 

Given  the  initial  conditions,  movement  between  the  states  will  continue  until  everyone  is  either 
vaccinated,  recovered  (or  immune)  or  dead.  The  initial  conditions  of  the  model  are  specified  by 
the  initial  number  of  people  infected,  and  the  initial  number  of  people  who  are  susceptible, 

U{to)-  The  total  number  of  people  in  the  simulation  is 

Nt^Uito)PE{to).  (3) 
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3.  SMALLPOX  MODEL 


In  the  smallpox  model,  the  six  primary  states,  shown  in  Table  1,  are  subcategorized  according 
to  whether  people  are  quarantined  or  isolated  (g),  unsuccessfully  vaccinated  designated  as 
having  a  traceable  meeting  with  an  infectious  person  (c),  receiving  treatment  or  will  die  {d). 
In  all,  there  are  twenty-one  states,  as  shown  in  Table  2. 


TABLE  2 


States  of  the  smallpox  model 


Disease  State 

Notation; 

Categorized  Disease  States 

Unexposed 

U{t) 

UvU) 

Uq{t) 

Uqv[t) 

Ucit) 

Ucv{t) 

Exposed 

E{t) 

Ev{t) 

Eq[t) 

Eqv[t) 

Ec[t) 

Ecv{t) 

Symptomatic 

Sit) 

Sd{t) 

Sqit) 

Sqd[t) 

Sqt{t) 

Sqtd[t) 

Vaccinated 

Vit) 

Recovered 

Rit) 

Dead 

D{t) 

3.1  PARAMETERS  OF  THE  SMALLPOX  MODEL 

The  parameters  of  the  model  can  be  grouped  into  four  categories:  (1)  disease  progression, 
(2)  prophylaxis  effiacy,  (3)  public  health  response,  and  (4)  resource  limits.  These  parameters  are 
presented  in  Tables  3-6. 

Table  3  lists  disease  progression  parameters  for  smallpox  up  to  and  including  the  infectious 
period.  Smallpox  has  a  relatively  long  incubation  period  followed  by  a  non-infectious  prodromal 
period,  marked  by  high  fever  and  flu-like  symptoms,  that  is  in  turn  followed  by  the  infectious 
period,  marked  by  the  characteristic  smallpox  rash.  Because  the  prodrome  is  assumed  to  be  non- 
infectious,  it  is  included  in  the  incubation  period.  The  final  step  to  recovery  is  progression  through 
a  non-infectious  scabbing  period,  which  is  not  modeled.  Thus,  R{t)  is  the  number  of  individuals 
who  have  survived  infection  with  immunity  but  who  may  not  be  fully  physically  recovered  yet. 

Currently,  there  is  no  treatment  to  improve  the  likelihood  of  survival  once  the  rash  has 
developed.  However,  we  have  included  vaxiables  related  to  treatment  (A^,  pt,  and  8t)  in  the  model 
so  that  the  model  may  be  applied  to  an  infectious  disease  that  does  have  effective  treatments  or  to 
model  the  outcome  of  a  smallpox  outbreak  if  an  acceptable  antiviral  drug  becomes  available. 

In  the  absence  of  any  medical  interventions,  whether  or  not  an  outbreak  becomes  an  epidemic 
depends  on  the  value  of  the  basic  reproductive  number,  i?o;  is  defined  as  the  expected  average 
number  of  secondary  infections  resulting  from  a  single  infectious  person  in  a  fully  susceptible  pop- 
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TABLE  3 


Smallpox  disease  progression  parameters 


Parameter 

Notation 

Value 

incubation  period  (days) 

14.6  [10] 

proportion  of  people  expected  to  die 

A 

0.30  [29] 

infectious  period  given  recovery  (days) 

8.5  [10] 

infectious  period  given  death  (days) 

8.5  [10] 

proportion  of  treated  people  expected  to  die 

0.30 

infectious  period  given  recovery  with  treatment  (days) 

Pt^ 

8.5 

infectious  period  given  death  with  treatment  (days) 

8.5 

proportion  of  contacts  infected 

4> 

0.20  [10] 

average  number  of  secondary  infections 

Rq 

10  [27] 

contacts  per  day  per  member  of  population 

see  equation  (4) 

general  population  (“herd”)  immunity 

e/! 

0.20  in  US  [27] 

environmental  infections  per  day 

7 

0 

Illation.  If  i?o  <  1,  then  the  outbreak  will  die  out.  However,  if  J?o  >  1  and  no  interventions  are 
imposed  to  force  Rq  <  1,  then  the  disease  will  continue  to  spread.  In  this  analysis,  the  contact  rate 
between  people,  6.  is  defined  as 


The  definition  of  (3  ensures  that  at  the  outset  of  the  outbreak  each  infectious  p)erson  infects  on 
average  Rq  people.  This  is  a  standard  way  to  define  the  contact  rate;  for  a  similar  definition 
see  [10]  in  which  contact  rate  is  defined  as  /3  =  Rq/ {(f>Nt6~^)  where  is  the  amount  of  time  an 
individual  is  infectious  [10]. 


TABLE  4 

Smallpox  prophylaxis  efficacy 

Parameter  Notation  Value 

prophylaxis  efficacy  when  uninfected  0.975  [10] 

prophylaxis  efficacy  when  infected  but  latent  te  0.30  [10] 


Table  4  lists  the  assumed  efficacy  of  the  smallpox  vaccine.  The  prophylaxis  efficacy  for  those 
already  infected  is  an  approximation  based  on  the  fact  that  the  vaccine  should  be  effective  at 
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preventing  or  reducing  the  severity  of  the  illness  if  administered  early  during  the  incubation  period, 
which  may  last  between  9  to  14  days  [25].  The  window  of  opportunity  for  effective  vaccination  is 
estimated  to  be  up  to  4  days  after  infection  [14]. 


TABLE  5 

Smallpox  public  health  response  parameters 


Parameter 

Notation 

Value 

time  elapsed  before  response  is  initiated  (days) 

K 

varied 

individuals  per  1  million  given  prophylaxis  per  day 

U 

100,000  [28] 

fraction  of  contacts  traced 

Be 

varied 

quarantine  period  (days) 

16.7  [10] 

average  time  delay  to  self-isolate  or  seek  treatment 

Bs 

0.95  [10] 

TABLE  6 

Resource  limits  for  a  smallpox  outbreak 


Parameter 

Notation 

daily  prophylaxis  available 

Mv 

total  prophylaxis  available 

^^'0  total 

daily  new  contact  quarantine  capacity 

M, 

daily  total  contact  quarantine  capacity 

total 

daily  treatment  capacity 

Mt 

Recently,  there  has  been  much  debate  about  the  proper  response  to  a  smallpox  outbreak. 
See  [6,11,18,20,35]  for  examples  of  mathematical  models  currently  being  used  to  compare  control 
policies  for  a  smallpox  outbreak.  In  this  paper,  we  consider  three  basic  responses:  (1)  contact 
tracing  and  ring  vaccination,  the  final  strategy  used  during  the  smallpox  global  eradication  cam- 
pedgn  [7,21];  (2)  mass  vaccination,  like  that  used  during  the  last  mass  vaccination  campaign  in  the 
US  [34]  (which  also  included  contact  tracing  [33]);  and  (3)  these  two  responses  in  combination.  In 
addition,  the  model  can  also  simulate  isolation  of  infectious  people,  which  would  also  be  another 
method  of  controlling  an  outbreak. 

In  the  event  of  a  smallpox  outbreak,  the  current  plan  based  on  CDC  guidelines  is  to  use  an 
‘'enhanced  version  of  the  surveillance-containment  strategy  that  was  successfully  employed  by  the 
World  Health  Organization  (WHO)  to  eradicate  smallpox  worldwide  in  the  sixties  and  seventies" 
[28]  and  to  resort  to  mass  vaccination  only  if  the  surveillance-containment  strategy  does  not  appear 
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to  be  effective.  In  the  case  of  mass  vaccination^  the  plan  is  to  vaccinate  100,000  people/day 
per  million  people  in  the  area  [28].  As  of  July  2003.  the  Massachusetts  Emergency  Management 
Association  plan  is  to  vaccinate  the  entire  state  of  Massachusetts  (6.4  million)  within  6  days, 
vaccinating  80%  of  the  population  in  the  first  3  days  and  the  remaining  20%  in  the  next  3  days  [23]. 
One  interesting  point  of  reference  is  the  last  mass  vaccination  campaign  that  occurred  in  the  U.S.  in 
1947,  which  began  after  a  man  with  undiagnosed  hemorrhagic  smallpox  infected  twelve  others  [7]. 
In  this  case,  6.35  million  people  in  New  York  City  were  vaccinated  between  April  4  and  May  2  [34] 
(although  there  is  some  debate  as  to  the  total  number  vaccinated  [33]). 

Yet  another  unknown  is  the  fraction  of  infected  people  who  would  be  able  to  be  traced  by 
the  public  health  system.  When  modeling  the  Kosovo  outbreaJc  of  1972  that  oc.curred  in  rural 
Yugoslavia,  Gani  and  Leach  atssumed  that  97.5%  of  infected  population  was  traceable  [10].  Given 
our  modern,  mobile  world  population,  this  percentage  may  be  much  lower. 


3.2  INFECTION  AND  CONTACT  TRACING 

In  this  section,  we  present  the  equations  that  govern  movement  between  unexposed  and 
exposed  states  via  ''meetings^  between  susceptible  and  infectious  people  as  well  as  environmental 
exposures  (7  /  0).  In  addition,  we  also  identify  whether  the  meetings  are  traceable  {6c  7^  0)  or  not. 
Some  equations  deal  only  with  identifying  those  individuals  who  have  had  a  traceable  meeting  and 
will  be  subsequently  contacted  and  either  vaccinated  or  quarantined. 

In  the  model,  contact  tracing  is  only  applied  to  those  who  have  not  yet  become  symptomatic, 
namely  people  in  states  Uc,  Ucv,  Ec  and  Ecv.  These  states  contain  people  that  have  had  traceable 
meetings  (as  designated  by  the  letter  c  in  their  state  names)  and  will  be  contacted  if  and  when 
contact  tracing  is  initiated.  Until  such  time,  people  in  Ec  and  Ecu  will  simply  progress  to  a 
symptomatic  state  as  if  they  were  in  E  or  Ev.  Similarly,  people  in  Uc  and  Ucv  will  be  just  as 
susceptible  to  infection  as  people  in  U  and  Uv. 

Below  ai’e  the  equations  governing  movement  out  of  states  tb  Uv,  Uc\  Ucv,  E  and  Ev.  A 
schematic  of  these  transitions  is  shown  in  Figure  1.  For  convenience,  let  the  number  of  people  who 
can  infect  others  (that  is,  the  number  of  infectious  people  not  isolated)  be 


7  =:  F  +  Sd, 


(^) 


and  let  the  number  of  unexposed  people  susceptible  to  environmental  exposure  be 


(j  =  U  ^Uq^  Uv  b  Uqv  +  Uc  b  Ucv. 


(6) 


The  first  four  equations  below  represent  traceable  meetings  with  infectious  people  while  the 


last  two  equations  represent  untraceable  meetings. 


fiU,Uc) 

1 

1 - 1 

u 

II 

(7) 

fiUv,Ucv) 

1 

o 

II 

(8) 

fiU,Ec) 

=  Occix^iu 

(9) 

f(Uv.Ecv) 

=  ec4>i3iuv 

(10) 

fiU,E) 

(11) 

f{Uv,Ev) 

Uv  Uv 

=  {l-9,)clxPIUv  +  mmi-^'r,—) 

(12) 

(Note  that  when  this  model  is  implemented,  care  must  be  taken  when  U  becomes  small;  this  issue 
can  be  addressed  by  adding  a  small  positive  number  to  V\  one  on  the  order  of  floating  point 
precision.) 

Equations  (13)  and  (14)  below,  which  are  similar  to  equations  (11)  and  (12),  represent  infec¬ 
tions  that  result  from  a  meeting  between  an  infectious  person  and  an  unexposed  person  who  has 
had  a  prior  meeting  with  an  infectious  person  that  was  traceable.  (Of  course,  this  prior,  traceable 
meeting  did  not  result  in  infection  as  they  axe  still  in  Uc  and  Ucik) 


f[Uc,Ec) 
f{Ucv,  Ecu) 


Uc  Uc 
<jypiUc  +  mm{^y,—) 

[/  Ar 

lo  Jtt  .  ,Ucv  Ucv, 

(pp  I  U cv  +  min(  7, 


U 


At 


(13) 

(14) 


The  equations  below  are  similar  to  equations  (7)  and  (8)  but  represent  traceable  meetings 
between  infectious  and  already  infected  people. 

f{E,Ec)  =  OcpiE  (15) 

f{Ev,Ecv)  =  OrpIEv  (16) 


The  equations  below  represent  infections  while  in  quarantine  due  to  environmental  exposure. 

f{Uq,Eq)  = 

f[UqxpEqv)  =  min(1^7,^^)  (18) 

If  7  =  0,  then  no  one  will  become  infected  while  in  quarantine. 

Expressions  involving  the  min(*,  •)  function,  which  represent  environmental  exposures,  also 
aj)pear  in  equations  (11)-(14),  but  not  in  the  equations  governing  movement  into  Uc,  Ucv,  Ec 
and  Ecu  because  movement  into  these  states  implies  that  the  source  of  infection  was  tixiceable. 
In  the  current  model,  we  assume  that  environmental  exposures  are  not  traceable.  As  a  side  note, 
environmental  exposures  may  be  traceable  if  a  large  number  of  people  can  be  traced  back  to  a 
single  place  where  the  exposure  occurred.  However,  this  is  currently  unlikely  in  a  time  frame 
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Figure  L  A  schematic  of  the  transitions  between  states  due  to  infection  and  contact  tracing.  Movement 
of  people  who  have  had  traceable  meetings  is  indicated  by  an  orange  arrow,  while  movement  of  those  who 
have  been  exposed  via  an  untraceable  meeting  is  indicated  by  a  brown  arrow.  Note  that  although  individuals 
in  lie  and  Uev  become  infected  via  an  untraceable  m,eeting  (hence,  the  brown  airows),  they  retain  their 
"'traceable^'  status  and  are  m,oved  into  Ec  and  Ecv.  Individuals  in  Uc  and  Ucu  have  had  traceable  contacts 
with  an  infectious  person  in  the  past:  however,  that  contact  did  not  lead  to  an  infection.  Finally,  infections 
due  solely  to  an  environmental  exposure  (^7  /  Oj  are  indicated  by  a  purple  arrow.  A  blue  circle  indicates 
that  those  in  this  state  have  been  vaccinated,  while  a  green  circle  indicates  tho,t  those  in  this  state  are  in 
quarantine. 
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on  the  order  of  a  week  or  so.  A  recent  example  of  identifying  a  common  source  of  exposure  is 
the  cluster  of  SARS  cases  at  the  Amoy  Gardens  Estate,  where  approximately  300  people  were 
infected  [30].  In  this  case,  transmission  appears  to  have  been  the  result  of  both  close  contact  and  a 
leaking  sewer  system.  It  was  contact  tracing  and  subsequent  analysis  that  led  to  the  identification 
of  a  possible  environmental  source.  Finally,  the  min(-,  •)  expression  does  not  always  yield  an  ODE. 
For  example,  when  {Uq/At)  <  y{Uq/U),  the  number  of  people  leaving  Uq  and  entering  Eq  is 
f{Uq,Eq)At  =  {Uq/At) At  —  Uq,  which  is  not  a  differential  equation. 

3.3  SMALLPOX  DISEASE  PROGRESSION 


The  following  equations  govern  movement  from  exposed  states  to  symptomatic  states: 


liE.S)  = 

(1  -  X)aE 

(19) 

f{E,Sd)  = 

XaE 

(20) 

fiEc,S)  - 

(1  —  X)aEc 

(21) 

f{Ec,Sd)  = 

Xa  Ec 

(22) 

f{Ev,S)  = 

(1  -  X)aEv 

(23) 

f{Ev,Sd)  = 

Xa  Ev 

(24) 

f{EciKS)  = 

(1  —  X)aEcv 

(25) 

II 

o 

Xa  Ecv 

(26) 

fiEq.Sq)  = 

(1  -  X)aEq 

(27) 

f{Eq,Sqd)  = 

XaEq 

(28) 

f{Eqv,Sq)  = 

(1  —  A) a  Eqv 

(29) 

f{Eqv,Sqd)  = 

Xa  Eqv. 

(30) 

These  equations  distinguish  between  the  fraction  of  those  infected  people  that  will 
and  of  those  that  will  die,  A. 

recover,  (1  —  A), 

The  following  equations  govern  movement  from  symptomatic  to  the  immune 

and  dead  states: 

f{S,R) 

=  pS 

(31) 

f{Sq,R) 

=  pSq 

(32) 

fiSqt.R) 

=  PtSqt 

(33) 

fiSd,D) 

=  SSd 

(34) 

f{Sqd,D) 

=  SSqd 

(35) 

fiSqtd,D) 

=  StSqtd. 

(36) 

As  mentioned  earlier,  the  infectious  period  is  followed  by  a  non- infectious  recovery  period, 
which  is  currently  not  modeled.  The  parameters  p  and  6  determine  the  amount  of  time  individuals 
are  in  S  and  Sd  and  allowed  to  spread  disease.  Thus,  R{t)  is  the  number  of  individuals  who  have 
survived  infection  with  immunity  but  who  have  not  likely  fully  recovered  yet.  See  Appendix  A  for 
equations  that  can  be  used  to  model  the  non-infectious  recovery  period. 
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3.4  OUTBREAK  CONTROL  POLICIES 


Whether  or  not  someone  is  vaccinated  or  quarantined  during  a  particular  time  step  (that  is, 
0?  where  j  is  a  state  holding  either  vaccinated  or  quarantined  people)  depends  on  whether 
or  not  there  is  room  in  quarantine  or  vaccine  available.  See  Appendix  B  for  the  scheme  used  to 
impose  the  resource  limits  defined  in  Table  6. 

3.4.1  Quarantine  and  vaccination  of  non-symptomatics 

In  this  section,  we  present  the  equations  governing  quarantine  and  vaccination  of  those  not 
yet  symptomatic.  Since  the  vaccine  efficacy  window  is  estimated  to  be  only  about  the  first  four 
days  after  exposure,  we  do  not  model  the  vaccination  of  those  already  symptomatic. 

Whether  people  are  vaccinated  or  quarantined  depends  on  whether  or  not  contact  tracing 
has  been  initiated  and  whether  quarantine  capacity  or  vaccination  supply  limits  have  been  reached. 
Quarantining  individuals  with  no  symptoms  can  only  occur  if  they  have  been  designated  as  traceable 
contacts  {6c  7^  0).  We  consider  three  possibilities:  (1)  both  vaccination  and  quarantine,  (2)  only 
vaccination  and  (3)  only  quarantine. 

Let  V  be  the  populations  in  the  model  that  may  request  and  receive  vaccination;  that  is, 

V^U^E.  (37) 

If  both  vaccination  and  quarantine  can  occur,  then  the  following  transition  rates  are  used: 


fiU,V) 

=  min( 

V  ' 

(38) 

f{E.V) 

=  Ce  min( 

^u,E) 

V 

(39) 

f{U..Uv) 

1  minf— ;y,  LO 

V 

(40) 

f{E,Ev) 

=  (1  -  e,.) 

E 

inin(  —  u,  E) 

■y 

(41) 

fiUc,V) 

=  eu  U c 

(42) 

f{Ec,V) 

—  Sc  Ec 

(43) 

f[Uc,Uqv) 

—  (1  ~ 

1  Uc 

(44) 

f{Ec,  Eqv) 

=  (1  -  ec) 

Ec 

(45) 

fiUcvMqv) 

=  Ucv 

(46) 

,f{Ecv,  Eqv) 

—  Ecu. 

(47) 

A  schematic  of  these  transitions  is  shown  in  Figure  2.  In  this  case,  when  people  move  to  quarantine, 
they  are  also  vaccinated  if  they  have  not  been  vaccinated  already.  Equations  (42)- (45)  model 
vaccination  as  a  result  of  contact  tracing. 
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Figure  2.  A  sckeinafAc  of  the  transitions  allowed  when  mass  vaccination  and  contact  tracing  are  instituted 
to  control  an  epidemic.  The  dashed  arrows  indicate  that  movement  will  only  be  allowed  if  there  is  sufficient 
capacity.  The  blue  airows  indicate  movement  due  to  the  mass  vaccination  protocol  The  pink  arrows  indicate 
movement  due  to  contact  tracing.  Successfully  vaccinated  individuals  resulting  from  contact  tracing  or  from 
the  rna.^s  vaccination  campaign  move  into  V.  (To  place  these  transitions  in  context,  we've  also  shown 
connections  between  the  states  representing  new  infections  and  identification  of  those  who  would  be  contacted 
if  contact  tracing  were  initiated.  As  in  Figure  U  brown  airows  indicate  movement  of  those  who  have  been 
exposed  via  an  untraceable  meeting,  while  movement  of  those  who  have  had  a  traceable  meeting  is  indicated 
by  an  orange  arrow.)  The  color  coding  of  the  circles  is  the  same  one  used  in  Figure  1. 
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Figure  3.  A  schematic  of  the  transitions  allowed  when  contact  tracing  is  followed  by  quarantine.  The  color 
coding  of  the  arrows  and  circles  is  the  same  as  that  in  Figure  2;  pink  arrows  indicate  transitions  due  to 
contact  tracing  and  quarantining  if  there  is  capacity. 


If  only  vaccination  can  occur,  then  equations  (52)-(43)  and  the  equations  below  are  used. 


f{lJcMcv)  = 

/ [Ec,  Ecv)  —  (1  —  ee)Ec. 


(48) 

(49) 


In  this  case,  equations  (42),  (43)  and  the  two  above  model  vaccination  as  a  result  of  contact  tracing. 
(In  the  experiments  to  follow,  we  assume  that  quarantine  space  is  unlimited.) 

If  only  quarantining  can  occur,  then  equations  (52)-(55),  (46),  (47)  and  the  equations  below 
are  used. 


RUc.Uq)  = 
f{Ec,Eq)  = 


A  schematic  of  these  transitions  is  shown  in  Figure  3. 


Uc 

Ec 


(50) 

(51) 


In  all  cases,  people  move  out  of  the  quarantine  after  an  average  of  a  ^  days: 


f{Uq,U) 

II 

(52) 

f{Uqv,Uv) 

=  aUqv 

(53) 

f{Eq,E) 

— 

(54) 

f{Eqv,Ev) 

=  oEqv, 

(55) 

where  a~^  should  be  greater  than  the  incubation  period,  so  that  those  who  are  infected  and  in 
quarantine  will  not  be  released  back  into  the  general  population  but  rather  will  become  sick  while 
in  quarantine  and  thus  be  held  there. 


3.4.2  Isolation  and  treatment  of  symptomatics 

Other  rates  (or  /(i,  j))  may  be  nonzero  due  to  isolation  and  treatment  of  symptomatic  people. 
Once  the  epidemic  control  policies  are  initiated,  we  assume  that  symptomatic  people  will  be  isolated 
with  an  average  delay  of  6~^  days  between  symptom  onset  and  isolation.  If  there  is  no  treatment 
available  and  symptomatic  people  can  only  be  isolated,  then  the  following  rates  are  used: 

fiS^Sq)  =  e,S  (56) 

J[Sd,Sqd)  -  e,Sd.  (57) 


If  we  are  modeling  an  infectious  disease  for  which  there  is  treatment  and  unfilled  capacity, 
then  the  following  rates  are  used: 


ns.sqt) 

=  OsS 

(58) 

f{Sq,Sqt) 

=  Sq 

(59) 

f{SE  Sqt) 

=  9,{l-~X)Sd 

(60) 

f{SE  Sqtd) 

=  OsXSd 

(61) 

f{Sqd,  Sqt) 

—  (1  -  A)  Sqd 

(62) 

f{Sqd,  Sqtd) 

=  XSqd, 

(63) 

where  A  =  A^/A  is  the  ineffectiveness  of  the  treatment.  Currently,  if  only  supportive  care  is  given 
to  smallpox  patients,  then  we  assume  that  Xt  =  A;  that  is,  mortality  remains  constant  despite 
treatment.  However,  antiviral  drugs  are  being  developed  [12,19].  Also,  we  assume  that  the  mean 
transition  time  to  enter  treatment  while  in  isolation  is  1  day  (implied  by  the  lack  of  the  factor  9s 
in  some  of  the  equations  above). 


3.5  MODEL  RESULTS 

3.5.1  Model  validation:  comparison  to  Kosovo  1972 

To  validate  our  model,  we  compare  the  output  of  the  model  (number  of  cases)  to  data  from 
the  Kosovo  outbreak  of  1972  [7,21],  for  which  the  global  eradication  strategy  was  implemented. 


15 


We  used  the  parameters  given  in  Tables  3-5,  except  those  shown  in  Table  7,  and  we  place  no  limits 
on  the  amount  of  vaccine  available  (daily  or  total)  or  quaxantine  space.  The  initial  conditions  are 
Eq  =  1,  Nt  =  2.2  million  people  [7].  In  addition,  we  set  Mf  —  0,  that  is,  no  simulation  of  treatment 
during  the  symptomatic  infectious  period. 


TABLE  7 


Parameters  for  the  1972  smallpox  outbreak  in  Kosovo 


Parameter 

Notation 

Value 

time  elapsed  before  response  is  initiated  (days) 

K. 

45.5 

number  of  individuals  given  prophylaxis  per  day 

U 

(l-e„)(iV,/45) 

general  population  (“herd")  immunity 

(h 

0.50  [10] 

fraction  of  contacts  traced 

Oc 

0.975  [10] 

daily  prophylaxis  available 

M, 

u 

total  prophylaxis  available 

total 

Nf  (unlimited) 

daily  new  contact  quarantine  capacity 

Mq 

Nf  (unlimited) 

daily  total  contact  quarantine  capacity 

^Qtotal 

Nf  (unlimited) 

daily  treatment  capacity 

Mt 

0 

On  March  16  vaccinations  in  Kosovo  began,  two  days  after  smallpox  was  suspected  in  4 
patients  [7],  45  days  after  the  index  case  is  estimated  to  have  become  infected  (k=  45.5  days).  By 
the  end  of  April,  after  about  45  days,  nearly  95%  of  the  population  of  Kosovo  {Nf  —  2.2  million)  was 
vaccinated,  implying  an  average  vaccination  rate  of  (A^^/45)  people/day  [7].  This  rate  corresponds 
to  a  rate  of  (1  —  e/j)(W^/45)  nomimmune  people/day,  given  the  initial  placement  of  ChNf  people  in 
V  to  simulate  ‘herd’  immunity  (Vq  =  ^^d  Uq  =  {1  —  e^)Nt  —  Eq). 

As  shown  in  Figure  4,  there  is  good  agreement  between  the  record('d  number  of  symptornatir 
individuals  and  the  number  predicted  by  the  model  for  Rq  =  10.  If  we  increase  Rq  to  10.2,  we  get 
a  better  fit  to  the  number  of  s^unptomatic  Cctses  between  50  and  70  days;  however,  doing  so  leads 
to  an  overestimation  of  the  cumulative  number  of  cases. 

Although  the  model  does  not  capture  the  typical  waves  of  cases  or  generations  of  illiK'ss.  w(' 
see  that  the  ODEs  do  capture  the  average  transition  time.  For  example,  we  vsee  that  the  index  case 
becomes  sick  after  approximately  15  days. 

As  shown  by  the  data,  the  incubation  period  of  smallpox  appears  to  follow  a  unimodal  dis¬ 
tribution  rather  than  an  exponential  one,  which  necessarily  underlies  the  ODE  that  governs  this 
transition.  This  is  one  reason  to  move  to  a  slightly  more  complicated  model  that  can  use  general 
probability  distribution  functions  (not  just  exponential  ones)  to  define  the  transitions  b(‘tw(^cn  the 
states. 
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Figure  4-  A  cornpanson  of  the  nuw,beT  of  cases  from  the  1972  Kosovo  outbreak  and  number  of  casts  predicted 
by  the  ODE  model  when  there  is  a  4^-5- day  delay  in  initiating  the  control  policy.  Depending  on  the  value,  of 
/?o  the  model  can  be  m.ade  to  fit  the  steep  increase  in  cases  (Rq  -  10.2)  or  the  cumulative  number  of  cases 
(R,  =  10.0). 


3.5.2  Hypothetical  contemporary  outbreak 

To  simulate  a  contemporary  outbreak,  we  used  the  parameters  given  in  Tables  3-5  except  those 
shown  in  Table  8,  We  consider  an  initial  release  of  smallpox  that  infects  100  people  (J5o  =  100) 
out  of  a  total  population  of  =  3  million  people,  which  is  approximately  the  number  of  people 
in  the  Boston  metropolitan  area,  and  we  place  no  limits  on  the  amount  of  vaccine  available  (daily 
or  total)  or  quarantine  space.  (See  Appendix  C  for  the  combinations  of  the  parameters  required  to 
model  the  public  health  responses.) 

If  the  U.S.  health  care  system  responds  as  the  one  in  Kosovo  did,  that  is,  during  the  second 
generation  of  secondary  cases,  then  the  model  estimates  that  there  will  be  roughly  106,600  cases 
and  over  31,800  deaths,  as  shown  by  the  red  curves  in  Figures  5  and  6.  The  second  generation  of 
cases  would  be  expected  to  surface  around  day  45  because  of  the  assumed  average  incubation  period 
of  14.6  days.  (The  first  cases  will  appear  around  day  15,  the  first  generation  of  secondary  cases 
around  day  30,  and  the  second  generation  of  secondary  cases  around  day  45.)  If  the  public  health 
system  is  able  to  diagnose  and  respond  after  the  first  few  cases,  but  before  the  first  generation  of 
cases,  then  the  number  of  cases  and  deaths  will  be  bounded  by  the  green  and  the  blue  curves.  On 
the  other  hand,  if  we  were  able  to  detect  and  then  initiate  a  response  once  the  first  case  appears, 
then  the  model  estimates  there  will  be  about  890  cases  and  roughly  270  deaths. 

Figure  7  shows  that  the  number  of  people  in  quarantine  peaks  soon  after  the  control  policies 
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TABLE  8 


Parameters  for  a  contemporary  smallpox  outbreak 


Parameter 

Notation 

Value 

time  elapsed  before  response  is  initiated  (days) 

n 

15,  30,  45 

number  of  individuals  given  prophylaxis  per  day 

V 

100,000  (iV</10®)  [28] 

fraction  of  contacts  traced 

Oc 

0.20 

daily  prophylaxis  available 

My 

total  prophylaxis  available 

My  total 

Nt  (unlimited) 

daily  new  contact  quarantine  capacity 

M, 

Ni  (unlimited) 

daily  total  contact  quarantine  capacity 

total 

Nt  (unlimited) 

daily  treatment  capacity 

Mt 

0 

are  implemented,  as  expected.  On  the  other  hand,  the  number  of  symptornatics  peaks  many  days 
(a  few  more  than  the  incubation  period)  after  the  policies  are  implemented,  as  shown  in  Figure  8. 


Figure  5.  The  cumulative  number  of  people  infected  during  the  contemporary  smallpox  outbreak,  S{t)  ~ 
S{t)  +  Sd{t)  -h  Sq{t)  +  Sqd{t)  +  R{t)  +  D{t). 


Figure  9  shows  the  number  of  deaths  for  each  of  the  response  delays  considered  and  for  each 
control  policy:  (1)  the  combined  policy  (contact  tracing  with  9c  =  0.20  and  mass  vaccination  at  a 
rate  of  100,000  people/day  per  1  million  people),  (2)  mass  vaccination  alone,  (3)  contact  tracing 
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alosie  with  a  very  optimistic  value  for  the  percentage  of  infected  people  found  [6c  =  0.90),  and  (4) 
contact  traicing  with  a  much  lower  value  for  the  percentage  of  infected  people  found  [6c  =  0.20). 

As  shown,  if  a  response  is  initiated  very  shortly  after  the  first  case  appears,  then  there  is 
not  much  difference  between  the  four  policies  considered,  but  a  mass  vaccination  campaign  would 
likely  result  in  more  vaccine-related  deaths  (which  is  currenlty  not  modeled).  The  number  of 
vaccine-related  deaths  is  expected  to  be  1  to  2  per  million  people  vaccinated  [22].  However,  rapidly 
identifying  and  tracing  the  initial  ca^es  infected  during  an  attack,  who  are  not  yet  symptomatic 
or  just  becoming  symptomatic,  might  be  difficult  without  biosensors.  As  delay  increases,  mass 
vaccination  and  the  combined  policy  are  equally  effective  at  controlling  the  epidemic. 

Finally,  it  should  be  noted  that  the  contact  tracing  rate  will  likely  be  limited:  this  is  not 
simulated  in  the  model.  In  addition,  random  effects  at  the  beginning  of  an  epidemic  are  not  modeled 
by  this  deterministic  model  [8]  nor  are  social  networks,  which  may  facilitate  contact  tracing. 


Figure  6.  The  cumulative  number  of  deaths  resulting  from  the  contemporary  outbreak,  D[t). 
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Figure  1  The  number  of  people  in  quarantine  (not  including  symptomatic  people)  during  the  coideniporary 
outbreak,  Q{t)  ™  Eq{t)  +  Uq{t)  -h  Eqv{t)  F  Uqv{t). 


Figure  8,  The  number  of  symptomatic  people  during  the  contemporary  outbreak,  S{f)  =  S{t)-\-Sd{t)-\-Sq{t)F 
Sqd{t). 


20 


delay  (days) 


Figure  9.  A  comparison  of  the  num.ber  of  deaths  given  three  different  sw,allpox  control  policies.  For  the 
policy  of  contact  tracing  alone,  we  considered  two  different  values  for  9c.  0.20  and  0.90.  For  each  policy,  the 
deaths  for  a  15-day  delay  are  under  LOOO;  they  are  269,  270,  502,  and  944  '^dh  respect  to  order  shown  m 
the  legend. 
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4.  A  MODEL  FOR  INHALATION  ANTHRAX 


Inhalation  anthrax  is  a  noncontagious  disease  that  we  have  modeled  using  three  stages  (a 
prodromal  period,  a  brief  ‘’recovery”  period,  and  then  fulminant  illness)  to  which  different  medical 
treatments  and  capacities  may  be  applied.  This  model  is  different  from  the  contagious  disease 
model  in  that  (1)  there  are  no  quarantine  or  contact  states  and  (2)  there  ai'e  multiple  stages  of 
the  disease.  In  addition,  for  smallpox,  individuals  are  protected  from  infection  by  vaccination; 
for  anthrax,  individuals  would  likely  be  given  a  60-day  supply  of  prophylactic  antibiotics,  such  as 
ciprofloxacin  [15]. 

In  this  model,  the  six  primary  population  states  are  categorized  according  to  whether  a 
person  is  unsuccessfully  given  prophylaxis  (u),  receiving  treatment  (t),  will  recover  (r),  or  will  die 
(d).  Table  9  lists  the  nineteen  states.  Figure  10  shows  a  complete  state  diagram  for  this  model, 
which  is  reflected  by  the  transition  rates  defined  in  sections  4.2  and  4.3. 


TABLE  9 


States  of  the  inhalation  anthrax  model 


Disease  States 

Notation: 

Categorized  Disease  States 

Unexposed 

U{t) 

Uv{t) 

Exposed 

E{t) 

Ev{t) 

Symptomatic  stage  1 

Slr{t) 

Sl{t) 

Slvr{t)  Slv2{t)  Sltr{t)  Slt2(t) 

Symptomatic  stage  2 

S2{t) 

S2tr(t) 

S2t3{t) 

Symptomatic  stage  3 

S3it) 

S3tr[t) 

S3tdit) 

Immune 

Vit) 

Recovered 

Rit) 

Dead 

D{t) 

Individuals  with  inhalation  anthrax  initially  experience  nonspecific  flu- like  symptoms  for  a 
period  of  days  and  then  sudden  fever,  labored  breathing,  and  shock  over  a  time  period  as  short  as 
a  few  hours,  followed  by  death.  However,  it  has  been  noted  that  there  often  appears  to  be  a  brief 
period  of  recovery  before  the  final  stage  [17].  In  this  model,  stage  1  represents  the  initial  nonspecific 
flu-like  stage;  stage  2,  the  brief  recovery  period  between  the  prodromal  and  fulminant  stages;  and 
stage  3,  the  fulminant  and  final  stage. 

When  creating  the  model,  we  assumed  that  prodromal  patients  will  receive  the  same  pro¬ 
phylaxis  as  that  given  to  patients  with  no  symptoms.  Thus,  unexposed,  exposed  and  prodromal 
people  draw  off  the  prophylaxis  supply.  To  prevent  so-called  double  counting  or  someone  receiving 
prophylaxis  before  and  after  they  become  symptomatic,  we  define  states  Slv  and  Slvr.  If  individ¬ 
uals  have  received  prophylaxis  when  they  are  latent  (infected  but  not  symptomatic),  but  it  was  not 
successful  (that  is,  they  are  in  Ev),  they  are  prevented  from  receiving  additional  prophylaxis  by 
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Figure  10.  A  schematic  of  the  transitions  between  states  of  the  anthrax  model.  The  dashed  arrows  indu  ate 
that  movement  will  only  occur  if  related  resource  lim/its  have  not  been  exceeded. 

moving  them  into  and  Slirr  rather  than  51  and  5lr.  See  Appendix  D  for  an  alternate  version 
of  the  model  in  which  the  prodromal  treatment  is  independent  of  prophylcrxis  supply  and  rate  of 
distribution. 

Currently,  the  model  does  not  consider  the  dependence  of  the  incubation  period  and  other 
parameters  on  the  dose  of  anthrax  received.  The  model  could  be  updated  to  include  a  number  of 
■‘levels’'  of  illness  corresponding  to  various  levels  of  dose  received,  thus  requiring  additional  disease 
parameters,  but  at  the  moment  there  is  relatively  little  data  from  animal  studies  to  support  such 
a  model. 


4.1  PARAMETERS  OF  THE  INHALATION  ANTHRAX  MODEL 


The  parameters  can  be  grouped  into  four  categories:  (1)  disease  progression,  (2)  prophylaxis 
efficacy,  (3)  public  health  response,  and  (4)  resource  limits,  as  shown  in  Tables  10-12.  Most  of 
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the  disease  progression  parameters,  shown  in  Table  10,  were  defined  based  on  data  from  the  most 
recent  inhalation  anthrax  cases  in  2001. 

We  chose  an  average  incubation  time  of  4.5  days  [17].  The  incubation  periods  that  could  be 
estimated  from  the  most  recent  cases  are  shown  in  Figure  11.  It  is  possible  that  longer  incubation 
periods  would  have  been  observed  in  2001  if  antibiotic  prophylaxis  had  not  been  distribued  to  the 
at-risk  population  [3,4].  Much  longer  incubation  periods  were  observed  during  the  1979  Sverdlovsk 
outbreak  [13,  26]  as  shown  in  Figure  11.  The  average  incubation  period  of  this  outbreak  was 
12.2  days,  which  is  longer  than  estimates  of  1  to  6  days  for  'Tigh-dose  exposures’"  [31].  Potential 
explanations  for  the  longer  incubation  times  are  delayed  germination  of  spores  [13,26],  low  dose 
exposure  [31],  and  reaerosolization  of  spores. 


TABLE  10 


Inhalation  anthrax  disease  progression  parameters 


Parameter 

Notation 

Value 

incubation  period  (days) 

a~^ 

1-7  [5] 

untreated  mortality 

A 

historically  0.89  [13] 

untreated  recovery  time  (days) 

25.1  [17,26] 

length  of  stage  1  (days) 

3.5  [17] 

length  of  stage  2  (days) 

1.5 

length  of  stage  3  (days) 

0.5 

treatment  inefficacy,  stage  1 

Aj 

0.45  [17] 

treatment  inefficacy,  stage  2 

A2 

0.95 

treatment  inefficacy,  stage  3 

A'3 

0.99 

treated  recovery  time,  stage  1  (days) 

Pi' 

22.3  [17] 

treated  recovery  time,  stage  2  (days) 

Pt2 

23.3 

treated  recovery  time,  stage  3  (days) 

Pt2 

24.3 

treated  length  of  stage  1  (days) 

3.5 

treated  length  of  stage  2  (days) 

1.5 

treated  length  of  stage  3  (days) 

0.5 

According  to  the  2001  data,  patients  took  approximately  3.5  days  on  average  to  go  to  their 
doctor  after  becoming  sick  [17].  Thus,  we  estimated  the  initial  stage  to  last  about  3.5  days,  which 
is  in  keeping  with  th<‘  CDC  fact  sheet,  which  states  1  to  5  days  [5]. 

To  define  the  remaining  disease  progression  parameters  for  untreated  individuals,  we  used  the 
data  for  the  two  patients  who  were  sent  home  without  receiving  antibiotics  and  who  later  returned 
to  the  emergency  room  with  fulminant  illness,  dying  after  1  and  3  days  after  their  initial  visits. 
For  these  patients,  an  average  of  5.5  days  passed  between  onset  and  death.  Thus,  we  estimate  the 
average  length  of  the  last  2  stages  to  be  2  days.  Since  stage  3  is  relatively  short,  we  define  the 
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Fiqure  IL  On  the  right  are  the  incubation  periods  of  the  2001  U.S,  outbreak  that  could  be  estimated  [17]. 
On  the  left  are  the  observed  incubation  periods  of  the  1979  Sverdlovsk  outbreak  [26]. 


duration  of  stage  2  to  be  1.5  days  and  the  duration  of  stage  3  to  be  0.5  days.  We  define  untreated 
recovery  time  to  be  25.1  =  21  [26]  +  4.1  [17]  days  for  stage  1.  For  stages  2  and  3,  we  increased 
this  recovery  time  by  1  and  2  days,  respectively.  We  assume  that  unsuccessful  treatment  does  not 
shorten  nor  lengthen  the  duration  of  the  symptomatic  disease  stages. 

Table  11  shows  the  assumed  prophylaxis  efficacy.  We  assume  an  efficacy  of  90%  when  pro¬ 
phylaxis  is  started  during  the  incubation  period.  Significant  protection  against  death  was  achieved 
with  prophylactic  antibiotics  for  Rhesus  monkeys  before  discontinuation  of  the  antibiotics  30  days 
post-exposure  [9].  However,  it  is  not  likely  that  there'  will  be  100%  compliance  in  the  completion 
of  a  60-day  course  of  antibiotics  [35].  Among  the  2,000  postal  workers  that  were  advised  to  take 
prophylactic  antibiotics  for  60  days,  it  is  estimated  that  only  about  40%?  of  them  completed  the  full 
course  of  antibiotics  [16]. 

The  response  parameters  considered  in  this  model  are  shown  in  Table  12.  These  parameters 
include  response  delay,  prophylaxis  distribution  rate,  and  average  time  delay  to  entc'r  treatment. 
In  addition,  capacity  limits  can  be  imposed  on  the  treatment  states.  These  limits  are  shown  in 
Table  13. 


TABLE  11 

Inhalation  anthrax  prophylaxis  efficacy 

Parameter  Notation  Value 

prophylaxis  efficacy  when  uninfected  Cu  0.90  [9,35] 

prophylaxis  efficacy  when  infected  0.90  [9,35] 
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TABLE  12 


Inhalation  anthrax  public  health  response  parameters 


Parameter 

Notation 

Value 

time  elapsed  before  response  is  initiated  (days) 

K 

varied 

individuals  given  prophylaxis  per  day 

U 

varied 

average  time  delay  to  enter  treatment  (days) 

07^ 

1 

TABLE  13 

Resource  limits  for  an  anthrax  outbreak 


Parameter 

Notation 

daily  prophylaxis  capacity 

My 

total  prophylaxis  capacity 

My  total 

daily  treatment  capacity,  stage  2 

Mt, 

daily  treatment  capacity,  stage  3 

Mt, 

4.2  INHALATION  ANTHRAX  DISEASE  PROGRESSION 

If  all  infections  occur  within  a  short  time  period,  for  example,  a  few  hours,  then  it  is  reasonable 
to  define  the  number  of  exposed  people  by  the  initial  condition,  E{to)  =  Eq,  rather  than  modeling 
their  inovenient  from  unexposed  to  exposed,  as  we  did  for  smallpox.  Since  we  are  considering  an 
aerosol  release  of  anthrax,  this  is  the  approach  we  took. 

We  begin  by  defining  the  transition  rates  from  exposed  to  symptomatic: 


/(E.51)  =  XaE  (64) 

fiE.Slr)  =  {l-X)aE  (65) 

f{Ev,Sl)  =  XaEv  (66) 

f{Ev,Slr)  =  {l-X)aEv.  (67) 

Once  symptomatic,  without  treatment,  the  disease  progresses  according  to  these  equations: 

f{Slr,R)  -  pSlr  (68) 

/(51,52)  =  (ii5l  (69) 

f(S2,S3)  =  62  S2  (70) 

/(53,D)  =  (53  53  (71) 
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With  treatment,  the  disecise  either  progresses  or  results  in  recovery  according  to  these  equations: 


f{Sltr,  R) 

=  pt^Sltr 

(72) 

f{Slt2,S2) 

=  5t^Slt2 

(73) 

f{S2tr'R) 

=  pt2  S2tr 

(74) 

/(52t3,53) 

=  6t2S2t3 

(75) 

fiSMr.R) 

=  ptsS3tr 

(76) 

f{S3td.,D) 

—  St^S3td 

(77) 

4.3  MEDICAL  INTERVENTIONS 


In  this  section,  we  define  the  transition  rates  related  to  prophylaxis  distribution  and  treatment. 
These  rates  are  non-zero  as  long  as  the  resource  limits  have  not  been  reached.  Let  V  be  the  states 
in  the  model  that  may  request  and  receive  prophylaxis;  that  is, 

K  =  +  51  +  51r.  (78) 


If  the  prophylaxis  supply  (daily  or  total)  has  not  been  reached,  then  the  following  equations  rates 
model  successful  and  unsuccessful  administration  of  prophylaxis  to  those  unexposed,  latent  and 
prodromal. 


ficr,Uv) 

fiE.V) 

fiE,Ev) 


■  ,E  E 

ep  rnm  —  /y,  (1  —  a)-— 

V 

E  E 

(1  -(e)  (1 


(79) 

(80) 
(81) 
(82) 


/{Sir,  Sltr) 

■  n 

=  mm(  ..  uAl  p)  ^  j 

E  At 

(83) 

f{Sl,Sltr) 

51 

=  (1  -  Ai)  min(-^ /A  (1  - 

5E 

(84) 

/(51,5U2) 

—  Ai  mm(-^  (1  -  Oi)  — ) 

V  At 

(85) 

Equations  (86)  and  (87)  define  the  rate  at  which  those  in  stage  2  and  3  enter  treatment  that  will 
ultimately  lead  to  recovery. 


f{S2,S2tr)  =  {l-X2)e,S2  (86) 

/(53,53tT)  =  (1-A3)6»s53  (87) 
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Equations  (88)  and  (89)  define  the  rate  at  which  thOvSe  in  stage  2  and  3  enter  treatment  that  will 
not  stop  the  disease  progression. 


f{S2,S2t3)  =  \2esS2  (88) 

f{S3,S3td)  =  XsOsSS  (89) 

4.4  RESULTS  FOR  A  LARGE-SCALE  AEROSOL  ANTHRAX  ATTACK 

The  model  was  used  to  investigate  the  effect  of  delaying  prophylaxis  and  treatment.  Of  course, 
the  effect  of  delay  can  be  heightened  or  lessened  by  decreasing  or  increasing  the  average  incubation 
time,  respectively.  We  chose  an  average  incubation  time  of  4.5  days,  based  on  the  most  recent 
cases  [17].  However,  as  mentioned,  the  Sverdlovsk  outbreak  did  have  a  longer  average  incubation 
time  [13,26]. 

For  the  large-scale  attack,  an  initial  number  of  exposed  people  was  determined  using  an  HPAC 
(Hazard  Prediction  and  Assessment  Capability^)  simulation  of  a  brief  line  release  of  anthrax  from 
a  boat  travelling  in  Boston  Harbor,  assuming  wind  from  the  south  east.  In  the  simulation,  1  kg  of 
ajithrax  was  released  in  16  minutes.  The  simulation  reported  that  approximately  510,000  people 
within  the  greater  Boston  area  would  be  exposed  (that  is,  receive  a  dose  greater  than  or  equal  to 
LCt2,  a  dose  lethal  for  greater  than  2%  of  those  exposed)  with  an  average  mortality  of  61%.  Thus, 
we  assume  that  510,000  people  are  infected  (£"0  =  510,000)  and  that  A  =  0.61. 


TABLE  14 


Parameters  specific  to  the  hypothetical  anthrax  attack 


Parameter 

Notation 

Value 

incubation  period  (days) 

4.5  [17] 

untreated  mortality 

A 

0.61 

daily  prophylaxis  capacity 

My 

100,000  -  750,000 

total  prophylaxis  capacity 

My  total 

Nt  (unlimited) 

daily  treatment  capacity,  stage  2 

Mt, 

0.50(16,600  [32]) 

daily  treatment  capacity,  stage  3 

Mt, 

600  [1] 

Next,  we  must  define  the  number  of  people  who  will  not  become  sick  but  will  demand  pro¬ 
phylaxis,  Uq.  If  we  define  this  number  to  be  small,  say  50,000,  this  corresponds  to  accurately 
identifying  those  exposed  and  preferentially  treating  them,  that  is,  having  very  accurate  knowledge 
of  the  shape  of  the  plume  as  a  function  of  time  and  those  people  in  the  plume  area  at  that  time. 
However,  if  there  is  no  way  to  determine  this  information,  then  with  so  many  people  becoming  sick, 
it  is  likely  that  the  entire  metropolitan  population  will  seek  prophylaxis,  and  the  total  number  of 

Mittp://www.dtra.mil/Toolbox/Directorates/td/programs/acec/hpac.cfm 
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people  in  the  simulation  who  will  demand  prophylaxis  should  be  reflective  of  the  number  of  people 
in  the  metropolitan  area  =  3,000,000,  Uq  =  Nt  —  Eq). 

For  this  scenario,  we  use  the  parameters  shown  in  Tables  10-12,  except  for  those  shown  in 
Table  14.  We  define  the  capacity  limit  for  those  in  stage  2  {Mt^)  to  be  a  percentage  of  the  number 
of  hospital  beds  in  Massachusetts,  which  is  estimated  to  be  16,600  [32];  we  assume  that  50%  of  these 
beds  can  be  made  available  to  patients  needing  advanced  care.  The  capacity  limit  for  stage  3 
is  estimated  by  the  number  of  ventilators  (portable  and  stationary)  in  the  push  packs  delivered  to 
New  York  City  on  September  11th  [1];  we  assume  that  existing  ventillators  in  the  city  will  already 
be  in  use.  It  should  be  noted  that  the  capacity  for  stage  2  and  stage  3  treatments  are  not  drivers 
of  the  attack  outcome  because  the  efficacy  of  these  treatments  are  so  low. 

Figures  12  and  13  show  the  total  number  of  deaths  and  symptomatic  people  as  a  function 
of  delay  in  public  health  response,  which  is  primarily  distribution  of  prophylaxis.  Each  pair  of 
colored  curves  corresponds  to  a  different  distribution  rate  {i/).  Figure  12  shows  that  as  treatment  is 
delayed,  the  number  of  symptomatic  people  approaches  the  number  of  those  exposed,  while  Figure 
13  shows  that  as  delay  increases  the  number  of  deaths  approaches  311,300  people  (or  XEo)- 

In  addition,  from  these  figures  we  can  see  that  the  number  of  cases  and  deaths  is  also  depen¬ 
dent  on  the  prophylaxis  demand.  During  the  incubation  period,  those  exposed  and  unexposed  will 
be  indistinguishable  from  each  other,  so  that  if  many  uninfected  people  demand  prophylaxis,  then 
those  who  are  infected  will  not  receive  prophylaxis  in  the  most  timely  fashion  possible.  Overall, 
the  sooner  medical  interventions  are  provided  and  the  higher  the  distribution  rate,  the  greater  the 
reduction  in  deaths. 

Assuming  a  higher  prophylaxis  efficacy,  for  example,  95%  due  to  higher  compliance  or  more 
effective  prophylaxis,  does  not  have  much  of  an  effect  on  the  number  of  cases  and  deaths  if  prophy¬ 
laxis  must  be  distributed  to  the  metropolitan  population.  Figures  14  and  15  show  the  predicted 
numbers  of  cases  and  deaths  assuming  95%  efficacy.  There  is  a  greater  difference  between  the  95% 
and  90%)  prophylaxis  efficacy  when  plume  mapping  can  be  used  to  target  prophylaxis  distribution 
(dashed  curves). 
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Figure  12.  The  number  of  symptomatic  people  as  a  function  of  delay  in  initiating  a  public  health  response, 
assuming  a  prophylaxis  efficacy  of  90%.  This  figure  shows  four  possible  distribution  rates  for  prophylaxis. 
Solid  lines  correspond  to  the  assumption  that  the  m^etropoUtan  region  will  seek  prophylaxis,  while  the  dashed 
lines  correspond  to  only  those  in  the  plume  plus  an  extra  50,000  people  dem.anding  prophylaxis. 


Figure  13.  The  number  of  deaths  as  a  function  of  delay  in  initiating  a  public  health  response,  assuming  a 
prophylaxis  efficacy  of  90%.  The  curves  are  colored-coded  according  to  prophylaxis  distribution  rate  as  in 
Figure  12. 
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550 


Figure  14-  The  number  of  symptomatic  people  as  a  function  of  delay  in  initiating  a  public  health  response, 
assuming  a  prophylaxis  efficacy  of  95%.  The  curtjes  are  colored-coded  according  to  prophylaxis  distribution 
rate  as  in  Figure  12. 


Figure  15.  The  number  of  deaths  as  a  function  of  delay  in  initiating  a  public  health  response,  assuming  a 
prophylaxis  efficacy  of  95%.  The  curves  are  colored- coded  according  to  prophylaxis  distribution  rate  as  in 
Figure  12. 
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5.  CONCLUSIONS 


The  disease  progression  models  can  shed  light  on  critical  windows  of  opportunity  for  public 
health  interventions  by  quantifying  the  benefit  of  early  response  and  high  treatment  rates.  In 
paxticulai',  if  antibiotic  prophylaxis  is  distributed  at  a  high  rate  within  the  first  day  following  a 
large  anthrax  attack,  nearly  all  those  exposed  can  be  saved  if  plume  localization  is  successful. 
The  results  indicate  plume  localization  is  most  beneficial  when  response  is  initiated  within  2  days 
post-attack. 

Similarly,  the  smallpox  model  indicates  that  early  initiation  of  a  mass  vaccination  campaign 
and  ciuarantine  can  significantly  reduce  the  number  of  deaths.  However,  the  current  smallpox  model 
may  underestimate  the  efTectiveness  of  contact  tracing  for  certain  outbreaks  since  it  does  not  model 
social  networks.  On  the  other  hand,  it  is  easy  to  imagine  a  threshold  number  of  index  cases  that 
would  overwhelm  the  contact  tracing  strategy  [28],  in  which  case  mass  vaccination  would  likely  be 
the  best  containment  strategy. 


33 


6.  FUTURE  WORK 


To  more  accurately  quantify  the  benefit  of  early  prophylcixis,  the  differential  equation-based 
model  should  be  updated  to  a  more  general  transition-based  model.  The  exponential  distribution 
underlies  the  transitions  of  the  ODE-based  model  and,  as  a  result,  leads  to  premature  transitions 
between  states,  in  paxticulax  the  transition  from  exposed  to  symptomatic.  For  example,  there 
were  no  cases  of  anthrax  2  days  prior  to  the  estimated  release  date  in  Sverdlovsk  [26],  Thus  the 
ODE  model  underestimates  the  benefit  of  early  intervention  since  treatment  is  less  effective  once 
the  individual  becomes  symptomatic.  An  updated  model  will  have  the  option  to  use  probability 
distribution  functions  defined  by  data  from  outbreaks  whenever  possible,  such  as  the  Sverdlovsk 
anthrax  outbreak  [26]. 
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APPENDIX  A 

MODELING  THE  NON-INFECTIOUS  RECOVERY  PERIOD 


Given  the  current  assumption  that  all  symptomatic  people  will  be  isolated  within  0^  ^  days, 
we  could  use  the  following  equations  to  estimate  the  number  of  physically  recovered  individuals, 
R{t). 


fiS,R) 

=  pS 

(A-1) 

fiSq,R) 

=  (P  +  Pn)Sq 

(A-2) 

f{Sqt,R) 

=  {Pt  +  Ptn)Sqt 

(A-3) 

In  this  case,  nearly  all  those  who  become  sick  move  into  R{t)  from  Sq{t)  and  Sqt{t)  with  average 
delays  of  p~^  ^  Pnt  respectively.  With  no  capacity  limits  imposed  on  isolation, 

only  a  fraction  of  symptomatic  people  would  move  to  R{t)  with  a  delay  of  p~^  days,  which  would 
essentially  represent  a  shortened  recovery  time.  However,  if  the  isolation  is  not  assumed,  then 
additional  states  and  equations  would  be  warranted.  In  particular,  a  state  to  hold  non-infectious 
recovering  people  who  are  never  isolated  would  be  required. 

In  the  case  of  limited  isolation  capacity  for  symptomatic  people,  one  approach  to  model  the 
number  of  fully  recovered  individuals  requires  one  additional  state  and  two  additional  parameters  (a 
stcite  to  hold  those  individuals  who  are  still  recovering  but  not  infectious  and  two  parameters  equal 
to  the  non-infectious  recovery  period  with  and  without  treatment).  We  only  need  one  additional 
state  for  two  reasons.  First,  we  can  assume  that  those  in  isolation  do  not  spread  disease  (given 
the  current  definition  of  /),  so  the  average  transition  time  between  Sq  and  R  can  be  increased  to 
the  infectious  period  plus  the  non-infectious  recovery  period,  and  similarly  for  Sqt  to  i?,  without 
increasing  disease  spread.  Second,  in  the  case  of  a  fatal  illness,  there  is  no  symptomatic  yet  non- 
infectious  period.  In  this  case,  the  equations  would  be 


/(5,5n) 

= 

(A-4) 

fiSn,R) 

=  Pr,Sn 

(A-5) 

f{Sq,R) 

=  {p  +  p.n)Sq 

(A-6) 

f{Sqt,  R) 

=  (pt  +  Ptn)Sqt, 

(A-7) 

where  Sn  is  the  number  of  recovery  people  who  are  no  longer  infectious,  Pn  is  the  non-infectious 
recovery  period,  and  pnt  is  the  non-infectious  recovery  period  with  treatment,  which  would  hopefully 
be  the  shorter  of  the  two  periods. 
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APPENDIX  B 

IMPLEMENTATION  OF  RESOURCE  LIMITS 


The  ODEs  presented  in  the  earlier  sections  were  solved  using  Euler’s  method  with  a  time 
step  =  0.1  day  prior  to  the  initiation  of  an  epidemic  control  policy.  After  the  initiation  of  this 
policy,  the  time  step  is  allowed  to  vary  so  that  the  capacities  for  quarantine  and  treatment  are  not 
exceeded.  To  do  so.  we  define 


(B-1) 


where 


AtM.^  -  {M,,-VdaUy{U))/{AV) 

~  Q  daily  /  i^Q) 

AtMt  =  i^t  -  Ty^atly{^^i))/{^T) 

~  ^^Qiotal  ~  QtotaliU))  I  {^Q)^ 


(B-2) 

(B-3) 

(B-4) 

(B-5) 

(B-6) 


where  AV,  AQ  and  AT  are  the  number  of  new  people  who  would  be  vaccinated,  moved  into 
quarantine  and  treated  (determined  by  evaluating  all  of  the  relevant  differential  equations  and 
cLSSuming  a  time  step  of  At  =  0.1);  Vdaiiy  is  a  running  total  of  all  prophylaxis  distributed  during  the 
current  day;  Qdatly  Tdaiiy  are  running  totals  of  people  in  quarantine  and  treatment  respectively; 
and  Vtotai  is  the  number  of  units  of  prophylaxis  given  out  so  far.  Note  care  must  be  taken,  when 
th(^  denominator  becomes  close  to  zero;  this  can  be  addressed  by  adding  a  tiny  number  to  the 
denominator. 


At  the  start  of  a  new  24-hour  period,  Vdaiiy  is  reset  to  zero.  However,  Qdmiy  and  Tdaiiy  are 
treated  differently;  people  must  recover  from  these  states  for  room  to  become  available. 

If  the  capacity  limits  for  any  of  the  interventions  is  not  met,  that  is, 


^^daily  i^'i) 

< 

(B-7) 

Q daily  {^1 ) 

< 

M, 

(B-8) 

'^daily  ) 

< 

Mt 

(B-9) 

^toial 

< 

^V,otai 

(B-10) 

Q total ) 

< 

total  ' 

(B-11) 

then  AE,  AQ  and  AT  are  calculated  as  follows: 

=  E  E  - E  (B-12) 

j  jitUqv  ]^Ea> 

where  j  =  [Uv  Uqv  Ucv  Ev  Eqv  Eev  V"]  is  a  vector  of  states  where  people  are  vaccinated  upon 
entrance.  Thus,  people  moving  from  Uv  into  Ev  and  Uqv  into  Eqv,  etc.,  are  not  counted  as  using 
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another  prophylcixis  unit.  (For  the  inhalation  anthrax  model,  j  =  \Uv  Ev  K].)  Similarly. 

where  j  =  [Uq  Uqv  Eq  Eqv  ]  is  a  vector  of  the  non-symptomatic,  quarantine  states;  thus,  people 
moving  from  Uq  into  Eq  and  Uqv  into  Eqv  are  not  recounted  as  entering  quarantine.  Also,  note 
that  quarantine  capacity  does  not  apply  to  those  who  are  already  symptomatic.  (For  the  inhalation 
anthrax  model,  quarantines  are  not  modeled.)  And  finally, 

AT  =  - ^  nSqtd.j),  (B-14) 

j  ji^St  J^Sqtd 

where  j  —  [St  Sid  Sqt  Sqtd  ]  is  a  vector  of  states  where  people  undergo  treatment  upon  entering 
these  states;  thus,  for  example,  people  moving  from  St  into  Sqt  and  Std  into  Sqtd  are  not  recounted 
as  entering  treatment. 
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APPENDIX  C 

RESPONSE  PARAMETER  VALUES 


In  the  next  two  paragraphs,  we  describe  the  parameter  values  required  to  initiate  the  public 
health  responses  we  are  considering  (using  the  Matlab  model):  (1)  contact  tracing  followed  by  quar¬ 
antine  and  vaccination,  (2)  mass  vaccination,  and  (3)  both  of  these  responses  in  combination.  To 
simulate  contact  tracing  followed  by  quarantine  and  vaccination  of  contacts  but  no  mass  vaccina¬ 
tion,  set  6c  G  (0, 1]  equal  to  the  proportion  of  contacts  expected  to  be  found,  set  ^ 

and  set  Mctotni  ^  >  0.  To  turn  off  mass  vaccination,  set  =  0.  If  />  —  0  and  ^  >  0^ 

then  only  those  contacted  will  be  vaccinated;  note,  or  the  mass  vaccination  rate  applies  only  to 
U  and  E. 

To  simulate  mass  vaccination  (of  people  in  U  and  E)  only,  define  the  vaccination  rate  >  0) 
and  set  ^  >  0.  (If  =  My  =  0,  no  one  will  receive  prophylcixis  even  if  /y  >  U.)  In 

addition,  set  dy  =  0  (so  that  people  in  U  are  not  moved  to  C/c,  a  state  to  which  mass  vaccination 
does  not  apply  and  similarly  for  E  and  Ec).  If  dc  —  0,  then  Uc{t)  =  0  and  Ec{t)  —  0  for  all  t  (and 
thus,  Ucv(t)  —  0  and  Ecv{t)  =  0  for  all  t);  this  ensures  that  all  those  unexposed  and  not  vaccinated 
will  remain  in  U  and  will  receive  prophylcixis  according  to  u  and  similarly  for  those  in  E. 
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APPENDIX  D 

ALTERNATE  ANTHRAX  MODEL 


In  section  4,  we  introduced  an  inhalation  anthrax  model  that  assumes  that  the  treatment 
for  prodromal  patients  is  the  same  unit  of  prophylaxis  distributed  during  the  mass  prophylaxis 
campaign.  In  this  model,  people  enter  the  prodromal  treatment  states  {Sltr  and  5ri^2)  at  a  rate 
of  i/.  We  could  assume  a  different  treatment  scheme  for  prodromal  patients,  either  the  number  of 
people  allowed  to  enter  the  state  can  be  rate- limited  using  a  different  rate,  or  the  number  of  people 
allowed  to  enter  can  be  based  on  a  capacity  limit  in  terms  of  people.  In  the  first  case,  the  capacity 
is  defined  by  the  number  of  people  treated  per  day  and  should  be  reset  daily.  In  the  second  case, 
space  in  the  state  becomes  available  as  people  recover. 

The  alternate  model  described  in  this  section  assumes  that  prophylaxis  and  the  prodromal 
treatment  ai'e  different,  so  states  Slvr  and  Slv2  axe  not  needed  to  ensure  that  people  will  not 
receive  prophylaxis  twice.  In  addition,  V  =  U  +  E.  Figure  D-1  shows  a  schematic  of  this  model. 
As  a  result,  i/  is  replaced  by  and  we  must  redefine  the  equations  governing  movement  into  the 
tr(!atmerit  states  Sltr  and  Slt2,  as  shown  below. 


f  {Sir,  Sltr)  =  dsSlr  (D-1) 

f  {SI,  Sltr)  =  (1-A;)0,51  (D-2) 

f{Sl,Slt2)  =  XiO.Sl  (D-3) 

Movement  into  Sltr  and  51^2  should  now  be  dictated  by  the  capacity  limit  Mt^,  and  capacity  will 
not  be  exceeded  by  appropriately  defining  At/  as 

min(At,  ,  AIm,  ,  ),  (D-4) 

where  AtMt  ~  AtM,.^)  and  A<,v/^ .  ,  and  At^i^  are  defined  taking  into 

consideration  the  respective  capacity  limits: 

{^III  -Tdail.yy{ti))/{ATi)  (D-5) 

AtM,^  =  {Mt,-T,iarly,iU))/iAr2)  (D-6) 

AtM,,  =  {Mt,-T,.aay,{k))/{AT3)  (D-7) 
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Figure  D~l.  A  schematic  showing  the  allowed  movement  between  states  of  the  alternate  anthrax  model  The 
dashed  arrows  indicate  that  movement  will  only  occur  if  related  resource  limits  have  not  been  exceeded. 
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